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Adaptive diffusion constrained total variation scheme with 
application to ‘cartoon + texture + edge’ image decomposition* * * § 
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Abstract 

We consider an image decomposition model involving a variational (minimization) problem and 
an evolutionary partial differential equation (PDE). We utilize a linear inhomogenuous diffusion 
constrained and weighted total variation (TV) scheme for image adaptive decomposition. An adaptive 
weight along with TV regularization splits a given image into three components representing the 
geometrical (cartoon), textural (small scale - microtextures), and edges (big scale - macrotextures). 
We study the wellposedness of the coupled variational-PDE scheme along with an efficient numerical 
scheme based on Chambolle’s dual minimization [14]. We provide extensive experimental results in 
cartoon-texture-edges decomposition, and denoising as well compare with other related variational, 
coupled anisotropic diffusion PDE based methods. 

Keywords: Image decomposition, total variation, linear diffusion, adaptive weights, multi-scale, 
denoising. 


1 Introduction 

Decomposing an image into meaningful components is an important and challenging inverse problem in 
image processing. Image denoising is a very well known example of image decomposition. In such a 
decomposition, the given image is assumed to be under the influence of noise, and the main purpose is to 
remove noise without destroying edges. This denoising task can be regarded as a decomposition of the 
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(a) Input (b) Cartoon (u) (c) Texture (v) (d) (Pseudo) Edges (w) 


Figure 1: MIT LabelMe decomposition result on an indoor image. Our scheme can obtain piecewise 
constant cartoon component along with texture and edges for better separation of basic shape elements 
present in a scene. Better viewed online and zoomed in. 

image into noise-free signal and noise part. There exist various methodologies for image restoration, where 
variational minimization and partial differential equation (PDE) are two of the most popular ones [3]. 

Another important example of image decomposition is based on its smooth and texture components 
using the total variation (TV) regularization method. This was first studied by Rudin et al [49] for image 
restoration. The TV regularization can be written as an unconstrained minimization, 

min E T y(u) = [ \Vu\ dx + — [ \u — f\ 2 dx. (1) 

u Jn 2a J n 

The parameter a > 0 balances the fidelity term with respect to the TV regularization. Let the given image 
be written as / = u + v, where the function u models well structured homogeneous regions (cartoon) and 
v defines oscillating patterns such as noise and texture. [34] established the scale connection property of 
TV for image decomposition, namely that the a parameter is related to the scale of objects in the image. 
In particular, Meyer proved that if / is a characteristic function and is sufficiently small with respect to a 
norm (| f\* < l/2a), then the minimization of the TV regularization (1) gives u = 0 and / = 0, which is 
counter-intuitive since one expects u = / and v = 0. Thus, Meyer proposed to use dual of the closure in 
the bounded variation (BV) space of the Schwartz class for non-trivial solutions, we refer to [34] for more 
details. This crucial fact has been exploited by [56] to obtain numerical approximations to the Meyer 
decomposition model, see also related early works [4,6]. 

Image denoising methods in general and variational PDE models in particular provide natural de¬ 
composition of images, see [3] for a review. Image smoothing and denoising literature is rich in terms 
of different methodologies utilized and current state of the art techniques [10,19, 43] have pushed the 
envelope of optimal improvements [18]. Next we provide a brief overview different approaches utilized in 
image decomposition. 
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1.1 Literature 


While the literature on cartoon and texture decomposition is extensive by now, three lines of inquiries 
are closely related to the work presented here. 

• Different function spaces for modeling the textures and discrete approximations: Fol¬ 
lowing Meyer’s seminal work [34], various authors have considered different function spaces to model 
textures accurately [5,6,24,27,29,31,32,40,54]. 

• L 1 fidelity based TV regularization models: In another related direction the fidelity term 
can be made L 1 and is proven to provide contrast preserving restorations. We refer to [16] for the 
geometric motivation and fine properties of L 1 fidelity term, see also [21]. Applications of Zd-TV 
models for cartoon plus texture decomposition are considered as well [2,22,27,37,58,59]. 

• Different Regularizers instead of TV: The well-known staircasing property of the TV regu¬ 
larization has been analyzed by many in the past [12,13, 38] and various modifications have also 
been studied. To avoid staircasing and blocky restoration results from the classical TV model 
there have been a variety of methods studied previously. Weighted and adaptive [45-47,51,52], 
nonlocal-TV [30], and higher order [7,17,28,41]. 

Out of other related decomposition models we mention multi-scale parameter based models [26, 53, 
54] which progressively capture texture. Extension to multichannel images in general, RGB images in 
particular, is also an important area of research [8]. Note that all of the above mentioned methods obtain 
cartoon and texture decomposition and in this work we obtain edges as part of our scheme. In Figure 1 we 
show an example decomposition for an indoor image taken from the MIT LabelMe [50] open annotation. 
Our advanced decomposition provides basic scene elements, since it is based on total variation (edge 
preserving) smoothing, the cartoon component reflects coarse shapes. This feature is useful in obtaining 
segmentations and annotations as the edges can provide guideline for separating salient objects [39]. 

1.2 Contributions 

In this paper, we propose a new image decomposition scheme which splits a given image into its ge¬ 
ometrical and textural part along with an edge capturing diffusion constraint. Following the success 
of weighted and adaptive TV models, our scheme is based on a weighted TV regularization where the 
edge-variable-based weight is computed in a data-adaptive way. The scheme is implemented using the 
splitting method of [9] along with dual minimization scheme for the weighted TV regularization [14]. As 
a by-product of the implementation via dual minimization we obtain an auxiliary variable which is akin 
to textural component of the given image. Thus, the scheme studied here provides a cartoon, texture, 
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edge (CTE) decomposition for digital images, see Figure 2 for some examples 1 . 

We consider the color image decomposition using the dual minimization based color TV model. Multi¬ 
scale decomposition following the recent work of [54] is also given. Moreover, we provide theoretical 
analysis of the model with a priori estimates and prove its solvability. Extensive experimental results 
on synthetic, textured real images are given. Further, illustrative comparison with other state-of-the- 
art models are undertaken and the results indicate superior performance of our scheme with respect to 
cartoon, texture separation as well as denoising with edge preservation. 

The rest of the paper is organized as follows. Section 2 introduces the adaptive TV regularization 
coupled with a diffusion PDE. Section 3 provides wellposedness results for our coupled model. In Section 4 
we provide the experimental results conducted on real and synthetic images and comparison with related 
schemes. Finally, Section 5 concludes the paper. 


2 Diffusion constrained regularization 


2.1 Weighted total variation minimization 

The total variation based regularization [49] given in Eqn. (1) is well-known in edge preserving image 
restoration, we rewrite it as follows, 


min El^ v 


(u) = f \Vu\dx + n f \u — fY 

Jn Jn 


dx 


( 2 ) 


where now g > 0 is the image fidelity parameter which is important in obtaining results in denoising and 
cartoon+texture decomposition. A related approach is to consider a weighted total variation, 


mm 


Eqtv( u ) = / g(x : u,\/u) \\/u\ dxg / \u — f\ 2 dx 

Jn Jn 


( 3 ) 


where g(x,u,Vu) represents the generalized weight function. For example, [9] have considered a convex 
regularization of the variational model 


lin < / g(x)\Vu\ dx + g / \u — f\dx 

u Un Jn 


using a fast minimization based on a dual formulation to get a partition in the geometric and texture 
information. Note that the image fidelity is changed to L 1 norm, we refer to [16] for more details. The 
convex regularization version is considered in [9] , 


min < [ g(x)\\7u\ dx -\—- [ (u + v — f) 2 dx + /i [ |u|dx>, 
u ' v Jn 2(9 J n J Q | 


( 4 ) 


x We use the words edges or pseudo-edges interchangeably in the text, since the edge component is computed with 
gradient maps without traditional edge refinement procedures such as maxima suppression [11]. 
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(b) Cartoon (u) 


(c) Texture v (d) (Pseudo) Edges (e) / — (u + v + w) 


Figure 2: The proposed coupled TV regularization with linear diffusion PDE model provides cartoon, 
texture and edge decomposition of images (stopping parameter e = 10 -4 ). Better viewed online and 
zoomed in. 



























(b) Cartoon component u with stooping parameters e = 10 4 , e = 10 5 , e = 10 6 and e = 10 7 using the proposed model 



(c) Cartoon component u with stooping parameters e = 10 4 ,e = 10 5 , e = 10 6 and e = 10 7 using the model in [9] 

Figure 3: (Color online) Our adaptive diffusion constrained total variation scheme (see Eqn. (9)) with 
constant /i, A (second row) provides better edge preserving image decomposition when compared to the 
traditional TV regularization model (Eqn. (4)) of [9] (third row) as stoping parameters increase from 
e = 10 -4 to e = 10 -7 . The proposed scheme keeps the structure without diffusing its boundary with the 
background. 

where the parameter 0 > 0 is chosen to be small so that / almost satisfies / ^ u + v, with the function 
u representing geometric information, i.e. the piecewise-smooth regions, and function v captures texture 
information lying in the given image. The function g is an edge indicator function that vanishes at object 
boundaries, for example, 

9 ^ X) 1 + /? |V/(x)| 2 ’ 

where / is the original image and /3 is an arbitrary positive constant. Thus, we see that TV based 
minimization models naturally lead to cartoon and texture decomposition of digital images. The image 
fidelity parameter /i can be made data adaptive to obtain texture preserving restorations, see [25]. 
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2.2 Description of the model 


In our work, we consider the following regularization model which was motivated by a coupled PDE 
modeling done in [48] for image restoration, 


mm 


(j)(x, u , | V^|) dx + /i L \u- f\dx ^, 


^ = \div(Vw) + (1 - A)(|Vw| - w). 


( 5 ) 

(6) 


The choice of regularizer depends on an application area and among a plethora of convex and non- 
convex functions available with the classical TV [49] and the non-local TV [30]. Motivated from the 
above discussions in Section 2.1, and success enjoyed by the weighted L X -TV regularization model in 
image denoising and segmentation, we use Zd-TV regularizer model as a prime example to illustrate our 
model here. The proposed CTE model thus consists of a minimization along with a non-homogeneous 
diffusion equation, 


min 

u 


g(w)\\7u\dx-\- J g{x) \u — f\ , 


( 7 ) 


dw 

dt 


= A (x)div(X7w) + (1 — A(x))(|Vu| — re), 


( 8 ) 


where g(w) = 1+ 1 w2 , or g(w) = exp(—w 2 ) (Perona-Malik type diffusion functions [42]), or g(w) = \w\ x , 

or g(w) = 7 1 (total variation diffusion function [49]). That is we solve adaptive data fidelity 

V e2 +M 

based weighted total variation minimization for the smooth part u using Eqn. (7) along with a linear 
non-homogenous diffusion constraint on w by solving Eqn. (8). Note that the balancing parameter A 
and image fidelity fi taking values in [0,1] are important in our experimental results. Adaptive ways of 
choosing these parameters are explained below in Section 2.3. Following [9] we use a splitting with an 
auxiliary variable v to obtain 


min 

u,v 


g(w)\Vu\dx + [ (u + v — f) 2 dx+ f g(x)\v\ dx >, 

Jn Jn 1 


(9) 


dw 

dt 


X(x)div(Vw) + (1 — A(x))(|Vr^| — w). 


Thus, the computed solution of these equations provides a representation (u,v,w), where the function u 
represents the geometric information, the function v captures the texture information, and the function 
w represents the edges lying in the given image. Figure 3 shows a comparison of our scheme (g = 1 and 
A = 0.5) and Bresson et al [9] scheme Eqn. (4) for a synthetic texture image which contains two different 
texture patterns. As can be seen, our scheme (Figure 3(b)) retains the cartoon edges better without 
diffusing the boundary and the shape is preserved in contrast to Bresson et al’s result (Figure 3(c)). 
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The above coupled system is solved in an alternating iterative way for all the variables (u,v,w) 
involved and Chambolle’s dual minimization scheme [14] is used for the weighted TV minimization step. 
We start with the initial conditions 

(u,v,w) | n=0 = (/, 0 , 1 ) 

and use the following steps to compute CTE components: 

1. Solving the linear diffusion PDE (8) for w with (u, v) fixed: 

w n+l =w n + \( x )Aw n + (1 - \(x)) (|V«| - W n )), (10) 

[dx) z 

where Sx is spatial discretization step (natural pixel grid), A is the standard finite difference dis¬ 
cretization for the Laplacian and St is the step size. 


2. Solving for the cartoon component u with (v,w) fixed: 

The minimization problem in u is given by (see Eqn. (9)), 

mini [ g(w)\Vu\ dx + [ (u + v — f) 2 dx\. (11) 

u Un 20 Jn J 

The solution of (11) is given by 


u = / — v — Odiv p, (12) 

where p = (pi,P 2 ) satisfies g(w)V(6 div p — (/ — v)) — \\7(0divp — (/ — v))\p = 0, which is solved 
using a fixed point method: p° = 0 and 

n+% _ P ra + 5tV(div(p n ) - (/ - v)/9) 

1 + ^R \V(div(p n )-(f-v)/6)\ 

3. Solving for the texture component v with (u, w) fixed: 

mm | ^ Jju + v - f) 2 dx + [i(x)\v\ drcj , (13) 

and the solution is found as 


v = 


/ -u- Ofi(x) 
< f — u + Ofi{x) 


0 

v 


if f — u > Ofjb(x), 
if f — u < —Q/jl(x), 
if \ f ~u\< 8[i(x). 


(14) 


Next we describe a data adaptive way for choosing the fidelity parameter /i using the cartoon component 
at a previous iteration u n . 


Remark 1 . We interchangeably use edges and pseudo-edges as the w component provides an edge like 
features from a given image. The definition of edges in a digital image depends on the context and many 
traditional definitions depend on the magnitude of gradients (i.e., \VI\), hence a solution of the PDE (8) 
provides a pseudo-edge map, see Figure 2(d). 







(a) Original image (b) /j, <(—adaptive, A ^—constant (c) fi constant, A ^—adaptive 


Figure 4: (Color Online) The proposed model (stopping parameter e = 10 -4 ) with adaptive fi\ (see 
Eqn. (17)) and constant A (= 1) provides better edge preservation in the cartoon component and captures 
small scale oscillations in the texture component against constant p (= 1) and adaptive A (using the 
definition of p i from Eqn. (17)). 


2.3 Adaptive fidelity parameter 

Here, we consider the data adaptive parameters selection strategies which can provide a balanced approach 
in obtaining better CTE decomposition results. For the image fidelity parameter p in Eqn. (7) we utilize 
a local histogram measure which was motivated from the image segmentation model of [36]. For a given 
gray-scale image I : 11 —> [0,L], let A f x , r be the local region centered at x with radius r. We compute 
the local histogram of the pixel x G ft and its corresponding cumulative distribution function 

\{z G A 4,r n ft I I(z) = y} I 


Px(v) = 

F x (y) = 


A 4,r n ft 

\{z G A r Xi r fl ft I I(z) < y}\ 
Af x r D ft 


for 0 < y < L, respectively. This allows us to define the following measurable function /i : ft —»> 
that for each x G ft, 


n(x) = m(x) = 


Jo FJjj)_dy 

m .ax Jq F x (y) dy’ 


(15) 

(16) 
such 

(17) 


allowing us to get a weight of how much nonhomogeneous intensity is present in a local region A f x , r of 
a given pixel x. This new feature of the image does not depend on the pixel properties instead provides 
regional properties, see [36] for more details. Thus, we see that the /i is chosen according to local histogram 
information and is computed in an image adaptive way using the cartoon u n in the iterative scheme (13). 
We compare our approach with two related adaptive functions: 

1. The adaptive formulation of [48] which uses a summation of cartoon components up-to iteration n. 


M2 (a:) = ^2 G Pi *w n (x) 


(18) 


i =0 


with pi = 1/i 2 and at n = 0 the \ 2 {x) = 0.05. 
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Figure 5: (Color Online) Comparison of different fi functions computed using the given input image, (a) 
Original image, (b) fii based on local histograms Eqn. (17). (c) /i 2 based on the work of [48] Eqn. (18). 
(d) /i 3 base on the work of [30] Eqn. (19). (e) Energy versus iteration for different adaptive pi functions 
based energy minimization scheme (7-8) with stooping parameter e = 10 -4 . 


2. Relative reduction rate based parameter proposed in local TV based scheme [30]. 


M3 0*0 = 


Gp * |V/(a;)| — G p * \L a * V/(a:)| 
G p * |V/(x)| 


(19) 


where L a is a low pass filter. Note that this adaptive parameter uses only the initial input image 
/ whereas the previous choices use u computed at a previous (Eqn. (17)) or every (Eqn. (18)) 
iteration. 


Remark 2. Figure 4 explains the choice of adaptiveness in our coupled model (7-8) for a synthetic image 
with different texture patterns. The first case with /a adaptive, X constant provides persisting cartoon 
component whereas the second case with /i constant, X adaptive (same local histogram based measure, 
Eqn. (17) is used for defining X(x)) blurs the boundaries in the final result. Thus, in what follows, we 
use only fi adaptive parameter to illustrate our decomposition results. 

Figure 5 shows a comparison of different adaptive pi functions for a synthetic texture image. We see 
that the local histogram based /ii captures the texture components from all the quadrants. Moreover, 
Figure 5(e) shows that the energy value decreases similarly for different /i functions as the iteration 
increases. Figure 6 shows the usage of different pi function when we apply our model (7-8) for the same 
synthetic image to obtain cartoon ( u ) + texture (v) + pseudo-edges (w) decomposition. Note that 
the texture image v is obtained by linearly transforming its range to [0, 255] for visualization purposes. 
Differences outside this range are saturated to 0 and 255 respectively. A similar transformation is applied 
to the edges (w) component as well. Next, we study the wellposedness of the model (7-8) using weak 
solutions concept and prove some a priori estimates and solvability of the proposed adaptive coupled 
model. 
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Figure 6: (Color Online) Different (i functions based CTE scheme (with A constant) results (stopping 
parameter e = 10 -4 ). Top: /i i (local histogram) based result. Middle: result. Bottom: result, (a) 

Cartoon u. (b) Texture v. (c) Pseudo-edges w. (d) v T w. Best viewed electronically, zoomed in. 


3 Wellposedness 


3.1 Preliminaries 


In the section, Q, is considered to be a bounded domain (i.e. an open set in R 2 ) possessing the cone 
property. We recall that this means that each point x G ft is a vertex of a finite cone C x contained in Cl, 
and all these cones C x are congruent [1]. Obviously, rectangular domains have this property. Fix also a 
time interval [0, T], T > 0. 

We study wellposedness of the weighted TV regularization, 


u(t , x) = u(x) : min 
u-.q^r 


g(w(t,x))\Vu(x)\dx + j fi(x)\u(x) - f(x)\dx L 


with the diffusion constraint, 

— = A p,a w(t,x) + (1 - X(x))(\Vu(t,x)\ - w(t,x)), 
w(t , x) = 0, x G dCl, 
w( 0, x) = F(x), 


( 20 ) 


( 21 ) 

( 22 ) 

(23) 


where p ^ 2, and f : Cl —y M, F : Cl —y [0, Too), A : Cl —y (0,1], fi : Cl — y (0, Too), g : [0, Too) — y (0, Too) 
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are prescribed functions. The operator A p ^\ is a weighted p-Laplacian: 

A p , x v = Adiv(|Vu| p_2 Vu) - (1 - \Vv\ p ~ 2 )\/v • VA. (24) 

In particular, for p = 2 we recover the linear diffusion case. In this section, for the sake of generality, we 
admit adaptive A and generic g. Note that w is non-negative by the maximum principle. We consider 
the Dirichlet boundary condition for w, but other boundary conditions can also be handled. 

We use the brief notations L q (q > 1), W™ (m G M), W™ ,q (m > 0) for the Lebesgue and Sobolev 
spaces on 12 with values in M or M 2 . Parentheses denote the bilinear form 

(u,v) = j u(x)v(z)dx. 
n 

The norm in L 2 is \\u\\ = yj(u, u). 

The symbols S(J\E), C(J\E), L\fJ\E) denote the spaces of Bochner measurable, continuous, 
Bochner integrable functions on an interval J C M with values in a Banach space E, respectively. 

Let M be the Banach space of finite Radon measures on f2. It is the dual of the space Co (Cl) (the 
space of continuous functions on 12 that vanish at 912, see e.g. [23]). 

Let BV be the space of functions of bounded variation on 12. For v G BV, and G C(i 2), 0 > 0, the 
weighted total variation of v is 

TV^v) = sup (n,div^)- (25) 

V » ecg °( n ):|^|<0 

In particular, the total variation of u is 


TV(v) = TVi(v). (26) 

Due to lower semicontinuity of suprema, for every non-negative <f G C{i 2) and a weakly-* converging 
sequence {n m } C BV, we have 

T^(r) < lim inf 7%(u m ). (27) 

m^+00 

A more refined argument of the same nature proves 

Lemma 1. For any cp G 5(0, T; C(£2)) ; cp > 0 for a.a. t G (0, T), and a weakly-* converging sequence 
{vm} C L q (0,T; BV), q > 1, one has 

2% (t) M*))< lim inf TV^ t) (v m (t)) (28) 

v ' m^+00 ' 

for a.a. t G (0, T). 

For v G BV , |Vu| will denote the corresponding total variation measure. The operator 

|V(-)I • BV M (29) 


12 



is bounded. We recall the duality relation 


TV^{y) — {\^v\,4>) MxCo (q.)- 


(30) 


The symbol C will stand for a generic positive constant that can take different values in different 
lines. We sometimes write C (..to specify that the constant depends on a certain parameter or value. 
We will use the embeddings 

BV C L q , q < 2, (31) 

W 2 C L q , q < + 00 , (32) 


and 


and the Poincare inequality 


Wl C C(Sl), p > 2, 


M C W~\ q < 2, 


IMIwi < C\\Vv\\ Lp , p > 1, v 

Embeddings (31)-(34) are compact (except for (31) with q = 2) 
We assume that A, VA and g are Lipschitz functions, 


dn 


= 0 . 


(33) 

(34) 

(35) 


A 0 = inf A(x) >0, g 0 = sup g(y) < +oo, 

xGO y>0 

and there exists a constant C g so that 

d(\ogg(y)) C g 

< -—for a.a. y > 0. 


dy 


1 + y 


The last condition means that g can have at most polynomial decay at infinity. 
We assume that /i G T^f}), and 


(36) 


0 < yi = ess inf y(x) < /i 2 = ess sup fi(x) < + 00 . 
x xeQ 


Finally, we assume that 

F G I/ 2 , 

and at least one of the following three conditions holds: 

P > 2, / G BV, 


p = 2,3q>l: feWj, 

1 

ff(y) 


p = 2, f E BV n Loo, 3c s > 0 : —— < c g (l + y), y> 0. 


(37) 

(38) 

(39) 

(40) 
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3.2 A priori estimates 


Before specifying the underlying function spaces and defining the notion of solution, let us derive a formal 
a priori estimate for problem (20)-(23). 

The Euler-Lagrange equation for (20) is 


t / / \ V u \ u ~ f r, du/dv 

— div giw)-—— + 11 - -- = 0, — . 

1 V J \\7u\J \u-f\ |Vtx| 


For each t G [0,T], multiplying (41) by w ^ w p , and integrating over H, we get 


= 0. 


(41) 


an 


Thus, 


Vu 

|V^ 


Ww(u - f) - 


f N \7u 

Wu g'(w)w 


|Vw| ’ g(w) 


Mu ~ f) 

g{w) 


Vw(u - f ) 


+ 


flW 

g(w) 


,\u- f\) =o. 


(42) 


Vu 


iv«r 


+ K|V«D- U^yv/I + ( ^,\u-f \) =0. (43) 


Multiplying (41) by and integrating over Q. we find 


Vw g\w ), 


Vw(u — /) ) + (1, |Vu|) — 


, |Vu| ’ g(w) 

Since the last term is non-negative, we conclude that 


Vu 

iv^i 


V/ 


g(w) 


g(w)' 


,\u- f\) =o. 


(44) 


TV(u) < TV(f) + C g \\Vw\\\\u - f\\. (45) 

Multiplying (41) by u — /, and integrating over 12, we derive 

(g{w), |Vu|) + \\n(u - f)\\ Ll < goTV(f) (46) 

Multiplying (41) by u|u|, and integrating over 12, we deduce 

2( fl HH,|V U |)+(|^i,H«l). (47) 

It is not difficult to obtain the following scalar inequality 

|a — b\ 3 < 2(a — b)a|a| + 2b 2 \a — fo|, a, b G M, (48) 

which enables to conclude from (47) that 

\Wg{ u - f) II 2 + 4(sMM, |v«|) < 2 II\/m/II 2 - (49) 

Hence, due to (31), 

||«-/||<C'(||/||)<C(||/||bv). (50) 
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By (45) and (50), 

TV(u)<C(l + \\Vw\\). 

Multiplying (21) by w , and integrating over £4, we get 


1 d\\w\\ 


+ (A, \Vw\ p ) + ((1 - A )w,w) = ((1 - A )w, |Vix|) - (wVA, \7w). 


2 dt 

Using Holder’s inequality and (36), we deduce from (43) that 

K |V«|) < (1 + C g )\\\7w\\ Lp \\u - f\\ Lp/f _ 1 + (w, | V/|). 

From (52), (53) and Young’s inequality we infer 

+ T fl Vu, llip - c (9Ao,p)\\u - /II v l v v ~\ + C(X,p)\\w\\ p L ^~\ + (w, | V/l). 

Provided (38) or (39), estimate (50), and embeddings (33) or (32), resp., imply 
d\\w\\ 2 


dt 


+ A 0 ||VHIL<C'(i + Hl p/p " 1 + Hkn- 


By (35),(37) and usual arguments, (55) yields 

IMUoofO,T;L 2 ) \\ W \\l p (0,T-,Wq' p ) — C- 

In case (40), we multiply (41) by and integrate over fl, arriving at 


(51) 


(52) 


(53) 


(54) 


(55) 


(56) 


Vu „ \ ( Vu q'(w)w „ , . .. 

uVw ) — [ — ——Vw I + (w, |Vu|) 


|Vu| : 


|Vu| ’ g(w) 




g(w)’ 


I u- fy g{w) 


= 0. (57) 


Then, since the penultimate term is non-negative, 


(w, |Vu|) < (1 + Cg)||Vtc|| ||u|| + fj, 2 Cg(\f\w, 1 + w) 


Now, (52) and (58) yield 


d\\w\\ 

dt 


+ Ao||V«;||L <C(l + ||tc|| 2 ), 


<C(1 + ||V W || + H| 2 ). (58) 


(59) 


which implies (56). In all the three cases, (46), (51) and (56) imply 


The operator 


IMIlooCO.TjL!) + \\u\\l p (0,T;BV) < C. 


A : BV M, 


(60) 
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A(u) = (l-A)|Vu|, 

(A(v), 4>)MxG 0 (n) = (|Vu|, (1 — A )4>)MxC 0 (n), 


is bounded, so 


||(1 - A)|Vu||| z , p ( 0) T;7W) < C. 

The weighted p-Laplacian operator 

is also bounded. Hence, (21),(34),(56) and (61) yield an estimate for the time derivative of w: 

q 


( 61 ) 


ll w IIl p (o,t;w,, -1 ) — q < 2, _ i >p. 


(62) 


3.3 Solvability 

Definition 1. Assume (37) and (38). A pair of functions ( u,w ) /rora the class 

u e Loo(0, T; Li) n L p (0,T; W), 

tc G Too(0,T; La) n L p (0, T; W 0 1 ’ p ) n W^(0, T; W^), 
is called a weak solution to problem (20)-(23) if 

TV g{w{t)) (u(t)) + || p(u(t) - f) || Ll < TV g(w(t)) {u(t)) + ||p(u(t) - f)\\ Ll , (63) 

for any u G 5(0, T; BV) and a.a. t G (0, T), 

w' — A p ^x + (1 — A )w = (1 — A)|V^| (64) 

in the space W~/ 1 p _ 1 for a.a. t G (0, T), and 

w(0) = F (65) 

™ w p ~/U- 

Remark 3. This definition is correct since all members of (64) belong to W~^_ x for a.a. t G (0, T) (cf. 
the end of Subsection 3.2), and w G W^(0, T; C C([0, T]; W r ^_ 1 ). 

Theorem 1. Assume (37) and (38). Then there exists a weak solution to (20)-(23). 

Proof. We can prove the existence of weak solutions via approximation of (20)-(23) by a more regular 
problem, and consequent passage to the limit (cf. [48,60]). Let (u m ,w m ) be a sequence of “approximate” 
solutions (possibly with “approximate” data f m and F m ). It essentially suffices to show that (63)-(65) 
is the limiting case of (20)-(23), i.e., that it is possible to pass to the limit in all the members. 
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Due to estimates (56), (60), (62), without loss of generality we may suppose that 


Um —^ u weakly — * in Loo(0,T; M), (66) 

u m -> u weakly — * in L p (0, T; BV), (67) 

w m w weakly — * in I/oo(0,T; L 2 ), (68) 

w m —>> w weakly in L p (0, T; W 0 1,p ), (69) 

w' weakly in L p (0,T; (70) 

Note that 

u G 7^(0, T; Lx) C 2^(0, T; M) fl L p { 0, T; BV). (71) 

By (33), (34) and the Aubin-Lions-Simon theorem [60], 


w m w strongly in L p (0, T; C(D)), (72) 

w m ^w strongly in C([0, T]; 1 ), (73) 

SO 

w m (0) ->■ w(0) in W~f p _ v (74) 

and we can pass to the limit in (65). 

Using the representation 

IMk = sup (<p,v), (75) 

veLoo, M Loo <i 

and lower semicontinuity of suprema, we can check that 


IIM«W - f)hi < 1™ inf \\fi(u m (t) — /)||i t 

m—)-+oo 

= lim inf \\/j,(u m (t) — f m )\\L 1 

m^+oo 

for a.a. t G (0,T). By Lemma 1, 

TV g{w{t)){u{t)) < m \un^miTV g{w{t ))(u m {t)). 

But 


(76) 


(77) 


k^s(wm(t)) TVg(w(t)) ( u m{t))\ 

< \\g(w m (t)) - g{w(t))\\ Loo TV(u m (t)) < C(g)\\w m (t) - w(t)\\ Lga TV{u m (t)), (78) 


so 


II TVg(w m ){u m ) - TVg( w )(u m ) ||l p/2 (0,T) < C\\w m - w\\ Lp ^0,T-,L oo )\\u m \\L p (0,T-,BV) 0. 


(79) 
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Therefore, without loss of generality, 


- TVg( w ( t ^(um(t)) -> 0 (80) 

for a.a. t G (0, T). Due to (76), (77), (80), we can pass to the limit in (63). 

On the other hand, (63) with u = u m , f = / m , w = w m , u = u gives 

TV g(w m (t))(u m (t)) + II n(u m (t) - fm) lUi < TVg{w m (t)){u(t)) + \\fJi(u(t) ~ fm)\\ Ll - (81) 

Similarly to (78)-(80), we can check that 

TV g(wm(t))(u(t)) - TV g ( w (t)){u{t)) -> 0. (82) 

From (76), (77), (80)-(82) we conclude that 


TVg(w(t)){ u m(t)) TV g ( w ( t )){u(t)) 


(83) 


for a.a. t G (0,T). 

Fix any non-negative function 0 G Co(0). Let 


K{t) 




Loo 


(84) 


and 

¥>(t) = K(t)g(w(t)) - 4>. (85) 

For a.a. t G (0, T), cp(t) is a non-negative continuous function on Q. By Lemma 1 and (83), we infer that 


TV<t>{u{t)) = k( t)TV g ( w{t )){u(t)) - TV v ( t )(u(t)) > lim sup(fc(t)TF 9(M(t)) (M m (t)) - TV v{t) (u m (t))) 

TYl t i"OO 

= lim sup TV^(u m (t)). (86) 

m—>-+00 


But, due to (28), 


2%(u(f))< lim inf TV^(u m (t)). 

ra—)>+oo 


Thus, 


7%(u(i)) = lim 7%(u m (f)), 

ra—)-+oo 


for every non-negative G Co(fi), which yields 


(87) 

( 88 ) 


|Vw m (f)| ->■ |Vu(f)| 


(89) 


weakly-* in M. for a.a. t G (0,T). Then (34) implies 


(1 - A)|Vw m (£)| —t (1 — A)|Vw(t)| 


(90) 
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strongly in W p / 1 p _ 1 for a.a. t G (0, T). Due to (67) and (34), 

||(1 — A)|Vrt m | \\ L T -w-. 1 ) — (91) 

5 J p/p— 1' 

By (90), (91) and [35, Proposition 2.8, Remark 2.10], 

(1 - A)|Vu m | ->■ (1 - A)| Vw|strongly in L q (0,T; Vq < p. (92) 

Rewrite (64) as 

w' + Aw = K(u, re), (93) 

where 


A(w) = - div(A|Vrc| p_2 Vre) + (1 - A )w, K(u , w) = -Vrc • VA + (1 - \)\\7u\. (94) 

It is easy to see that the operator A : Wq’ p —>> W p } p _ x is monotone, coercive and hemi-continuous 
(cf. [33]). By (72) and (92), 

K{um , Wm) -► K(u,w) strongly in L p / p _i(0, T; W^ 1 ^). 

Hence, we can successfully pass to the limit in (93) via Minty-Browder monotonicity technique (cf. 
[33]). □ 

Definition 2. Assume (37) and (39) or (40). A pair of functions (u,w) from the class 

u G Loo(0, T; Li) D L 2 (0, T; HR), (95) 

«;GL oo (0,r;L 2 )nL 2 (0,r;W 0 1 ’ 2 )nT^ 2 1 (0,T;W- 1 ), \fq < 2, (96) 

zs called a pseudosolution to problem (20)-(23) if there is a sequence (iz m ,rc m ,p m ) srzc/z that each pair 
(u m ,w m ) is a weak solution to (20)-(23) with p = Pm, 

u rn -G weakly — * in 1/^(0, T; At), 
u m —)> u weakly — * in 1/2(0, T; HR), 

Wm w weakly — * in 1/^(0, T; H 2 ), 

—>> re weakly in H 2 (0,T; IRq 1 ’ 2 ), 
re m w strongly in L 2 (0, T; L q ), Vg < Too, 
w m -> W strongly in C([0, T]; W^ 1 ), Vg < 2, 
re^ —>> w' weakly in H 2 (0,T; HR -1 ), Vg < 2, 

Pm —>• 2. 

Theorem 2. Assume (37) and (39) or (40). T/ien £/mre exists a pseudo solution to (20)-(23). 

Proof. The proof is based on estimates (56), (60), (62) and the proof of Theorem 1. □ 
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4 Experimental Results 


4.1 Implementation details 

The proposed scheme is implemented using the dual minimization [14] for the weighted TV (Eqn. (7)) 
and explicit Euler finite difference scheme for the non-homogenous linear diffusion (Eqn. (8)). The edge 
indicator function g(w) = 1/(1 + w 2 ) is used for all the results reported here. We obtained similar results 
for other g functions. The adaptive fi i based results are reported here unless otherwise stated explicitly 
and fi 2 provided similar results whereas provided slightly blurred cartoon components, see Section 2.3 
for details. The parameters Sx = 1, St = 1/8 and 6 = 10 -2 are fixed, and the best results according to 
the max(\u nJrl — ^ n |,|f n+1 — v n \< e) are shown. By constant A and g in the results we mean they are 
taken as constant value 0.5 and 1 respectively. The implementation of the proposed scheme is done for 
the constant choice (g e-constant, A e-constant), the adaptive choice (/i v-adaptive, g e-constant) and 
the multiscale case (/i multiscale, A constant). 

The algorithm is visualized in MATLAB 7.8(R2009a) on a 64-bit Windows 7 laptop with 3GB RAM, 
2.20GHz CPU. It takes on average <10 sec for 50 iterations for 3 channels image of size 481 x 321. 
Implementation is done over the following databases: Brodatz texture collection 2 , Mosaic art images 3 , 
Kodak Lossless True Color Image Suite 4 , Color Test Images Database 5 , USC-SIPI Image Database 6 and 
Simulated Brain Database 7 , Berkeley segmentation dataset of 500 8 . 

4.2 Image decomposition results 

4.2.1 Gray-scale images 

We first show decomposition results of [9] with our model in Figure 7 for a synthetic image which 
consists of two different texture regions. Comparing the cartoon - texture decomposition of our scheme 
(Figure 7(b)) with the results of Bresson et al (Figure 7(c)), we see that they behave different visually. 
For example, the shape of the diamond at the center is preserved well in our scheme whereas the [9] 
scheme blurs it in the final result. Figure 7(c) shows the energy value against number of iterations for 
the same synthetic image, which indicates that our adaptive CTE scheme decreases the energy values 
comparable to [9] model. More grayscale image decomposition results are given in Figure 2. We see 
that the cartoon component obtained are piecewise constant segments indicating the potential for image 

2 http: //mult ibandtexture. recherche, usherbrooke.ca 

3 htt p: / / www. cse .cuhk.edu.hk/leojia/ pro j ect s / text uresep / 

4 http://rOk.us/graphics/kodak 

5 http: / / www.hlevkin.com/TestImages / classic, htm 

6 http://sipi.usc.edu 

7 http://brainweb.bic.mni.mcgill.ca/brainweb/ 

8 http: / / www.eecs.berkeley.edu/Research/Projects/CS/vision/bsds / 
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(a) Original 


(b) u & v (Our) 


(c) u & v ( [9]) 



(d) Energy Vs Iterations 


Figure 7: (Color online) Comparison of our constant /q A proposed scheme (second column) with Bresson 
et al [9] (third column), shows that our scheme preserves large-scale textures and shape boundaries 
(stopping parameter e = 10 -4 ). Energy value comparison between our scheme CTE and Bresson et al. [9] 
shows similar convergence property. Best viewed electronically, zoomed in. 

segmentation [44]. The texture and edges component are complementary and it is clear that edges are 
based on the cartoon subregions, see for example, Figure 2(d) top row. 


4.2.2 Color images 

We further provide image decomposition for color images [8,22] by using vectorial TV version in Eqn. (7) 
following [8]. We consider the following vectorial TV with non-homogenous diffusion equation: 

min 

u=(ui,u 2 ,u 3 ) 


9(wi) 


\ 


3 3 » 

J2\Vui\dx + J2 / M( 

i =1 i= 1 n 


x) \ui —fi\dx 


dun 

dt 


Xi(x)div(Vwi) + (1 - Aj(a;))(|Vwj| - iUj), 


were each scalar function Ui : Q —► R, l<i<3 represent one component of the RGB color system. Note 
that following our alternating iterative scheme given by equations (10), (12) and (14) the new solutions 
(' Ui,Vi,Wi ) are given by: 

Ui = fi-Vi- Odiv (97) 
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(a) Original Image 


(b) Our constant p,, A based scheme 


(c) Bresson et al [9] 




(a) Original Image 


(b) Our adaptive pu, constant A based scheme 


(c) Bresson et al [9] 


Figure 8: (Color online) Our diffusion constrained total variation scheme (stopping parameter e = 10 -4 ) 
provides better edge preserving cartoon component u when compared to the traditional TV regularization 
model [9]. Even with constant /i, A the proposed scheme provides better results (see top row (b)). The 
crop regions highlight that the proposed scheme provides better preservation of large scale textures 
compared to [9] model. Best viewed electronically, zoomed in. 


where = (pi i: Pi 2 ) satisfies g(wi)\7{0 div — (fi — Vi)) — |V(#divp^ — ( fi — t^))|p^ = 0, which is solved 
using a fixed point method: p^ = 0 and 


P? +1 = 


Pi + 5t\7 (div(p ™) — (fi — Vi)/0) 


1 + 


St 


a(wi)\ jri 


ElV(div(P ?)-(/i-«i)/0)l s 


fi-Ui- 9/i(x) 


fi-Ui + 6n(x) 


Vi = 


0 
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if fi-Ui> 0 /j,(x), 
if fi-Ui < —9/j,(x). 
if | fi - Ui\ < 0n(x). 


(98) 

































Figure 9: (Color online) Effect of constant fi in color image decomposition using our coupled scheme 
on the cartoon (u) component. First row: Original color RGB images. Cartoon component result for 
constant /i, A with stopping parameter, Second row: e = 10 -4 , and Third row: e = 10 -6 . As can be seen 
decreasing the stopping parameter removes more texture details and provides piecewise constant cartoon 
image. Last row: Shows the proposed scheme results with adaptive //i, see Eqn. (17) and constant A 
(stopping parameter e = 10 -4 ) . Best viewed electronically, zoomed in. 
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Figure 10: (Color online) Corresponding texture components v. Arrangement is as in Figure 9. Best 
viewed electronically, zoomed in. 



Figure 11: (Color online) Corresponding edge functions w. Arrangement is as in Figure 9. Best viewed 
electronically, zoomed in. 
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< +1 = K + ^2 (A(*)A< + (1 - A(x)) (|V Ui | - O). (99) 

Figure 8 shows the cartoon components of our CTE scheme with constant and adaptive /ii against 
the traditional TV based scheme of [9] for two standard RGB test images. As can be seen, our scheme 
obtains better edge preserving cartoon (u) components. The close-up shots indicate that our scheme also 
avoids blurring of edges, see for example Barbara face. See also Figure 3 where the proposed scheme 
with adaptive /ii provides better shape preservation as stopping parameters increase from e = 10 -4 to 
e = 10 -7 . Next, Figures 9-11 shows decomposition for a variety of RGB images for two different iteration 
values in our proposed CTE scheme with constant fi against adaptive /ii based results. As can be seen 
in Figure 9, increasing the number of iterations removes more texture details and provides piecewise 
constant (smoother) cartoon images. Our adaptive /ai based scheme results (last row) on the other hand 
keep most of the salient edges. This can be seen further in Figure 10 (last row) where the adaptive fi i 
based scheme captures only small scale oscillations corresponding to texture components whereas the 
constant (i based results remove strong edges as well. Figure 11 show the corresponding edge functions, 
and it can be seen that adaptive scheme has more information than in constant parameter case. Thus, 
we conclude that, using adaptive /i provides an image adaptive way of obtaining edge preserving cartoon 
components without sacrificing overall decomposition properties. 

Figure 12 shows the cartoon component of two different color images as well as the edge functions 
given by the proposed model, and Canny edge maps [11] computed in all the three channels 9 . As observed 
the function w captures edges based on the cartoon subregions i.e. large intensity gradients while the 
Canny edge maps captures large and small intensity gradients. 

4.2.3 Comparison with other decomposition methods 

Next we use the grayscale digitize images from USC-SIPI Image Database to compare the performance 
of the CTE model together with the adaptive and multi-scale case with the following schemes: TV- 
47 -1 decomposition model (OSV, using finite differences) [40], T \-G v decomposition model (VC, finite 
differences channel-wise) [55], T V-Garbor decomposition model (AGO, Chambolle’s projection algorithm 
) [6], BV 2 — L 2 model (BP, which uses an algorithm close to the one set by dual minimization of [14] 
and [55]) [7], TV — L 1 model (CP, a first-order primal-dual algorithm) [15], 17-TV scheme for color 
images (DHC, using a unified Moreau-Yosida based primal-dual approach) [20], and the NL- Means 
algorithm (BCM, a improvement of the classical and state-of-the-art methods in image denosing) [10] 
which computes the local average color of the most resembling pixels using non-local means. 

Numerical results are given in Table 1 for sixteen standard natural grayscale images of size 256 x 

9 Computed using the MATLAB’s default command edge(u,‘canny 5 ). 
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(a) Input (/) and u components using the proposed model 



(b) (Pseudo) Edges (w) (c) Canny Edges 


Figure 12: (Color online) Pseudo edge maps & Canny Edges. First Row: Input (/) and u components 
using the proposed model with (/i adaptive, A ^—constant) and stopping parameter e = 10 -4 . Second 
and Third Row: (Pseudo) Edges (w) (First Column) and Canny Edges (Second Column). 
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USC-SIPI database with size f = 256 x 256 (Noisy PSNR = 22.11) and J = 512 x 512 (Noisy PSNR = 22.09). Each row indicates PSNR/MSSIM 
values for different test images. The diffusion with multiscale {/jl multiscale, A constant), adaptive choice (/i adaptive, A constant) 
and constant CTE scheme (/i ^—constant, A ^—constant) are given as MCTE, ACTE and CCTE (last three columns) respectively (stopping 
parameter e = 10 -4 ). Compared methods are 0SV [40], VC [55], AGO [6], BP [7], CP [15], DHC [20], BCM [10]. Best results are indicated in 
boldface and the second best is underlined. 





































(a) Barbara image (b) Canny edges (c) Noisy Barbara (d) Canny edges (e) Input noise 




(e) u, v, w, f — u , and / — (u + v + w) components for the proposed scheme. 


Figure 13: Decomposition of noise inputs, (a) House , F — 16 and Barbara images corrupted by Gaussian 
noise images (standard deviation a = 30). (b)-(d) u & f—u components provided by the VC [55], OSV [40] 
and AGO [6] schemes, respectively, (e) u,v,w, f — u and / — (u + v + w) components using the proposed 
CTE scheme with (CCTE - /a e-constant, A e-constant) for the Third row, with (ACTE - fi adaptive, 
A e-constant) for the Fourth row and (MCTE - fi multiscale, A <—) for the Fifth row. The corresponding 
PSNR (dB) and MSSIM measures for each model are displayed in Table 1. 
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256 and 512 x 512. As we can see the proposed scheme based on multi-scale lambda (/i <— multiscale, 
A ^—constant) and the adaptive choice (/i multiscale, A constant) performs well for a variety of images. 
Even with the adaptive choice (ACTE, (i adaptive, A e-constant) outperforms previous variational 
models bases if -1 , G p and Garbor norms and in some cases the L T -TV and the NL-Means schemes. 
For textured images such as Mandrill and Barbara we perform all variational models except the NL- 
means algorithm. We remark that the proposed CTE model does not aim to give state-of-the-art results 
for image denoising, and instead concentrates on demonstrating how our decomposition model can be 
harnessed for noise removal and edge detention. Figure 13 shows decomposition for Barbara into different 
components for some of the schemes in Table 1. We note that for other schemes we are able to discompose 
the image in u (smooth) and / — u (random noise). In our case, we are able to obtain edge variable w 
component. We also notice that for our model the random noise component f — u is directly given by 
the v component. 


4.3 Image denoising results 

4.3.1 Error metrics 


Currently there are no quantitative ways to evaluate different decomposition algorithms. In particular 
which smooth, texture and edge separation model are better is an open research question in image quality 
assessment. The proposed decomposition provides piecewise cartoon component which is obtained using 
a weighted TV regularization in an edge preserving way, see Figure 9 (last row). Hence, as a byproduct we 
obtain image denoising, with u the denoised image and v + w the ‘noise’ part. To compare the schemes 
quantitatively for the purpose of denoising, we utilize two commonly used error metrics in the image 
denoising literature, one is the classical peak signal to noise ratio (PSNR) [3], and the other one is the 
mean structural similarity measure (MSSIM) [57]. 


1. PSNR is given in decibels (dB). A difference of 0.5 dB can be identified visually. Higher PSNR 
value indicates optimum denoising capability. 

PSNR(m) := 20 * log 10 dB, (100) 

where MSE = (mn) -1 — uq), m x n denotes the image size, u max denotes the maximum 

value, for example in 8-bit images u max = 255. 


2. MSSIM index is in the range [0,1]. The MSSIM value near one implies the optimal denoising 
capability of the scheme [57] and is mean value of the SSIM metric. The SSIM is calculated 
between two windows uq and Co >2 of common size N x V, 


SSIM(cji,cj 2 ) 


(2/i UJl du 2 + c l)(2cr + c 2 ) 

(/4i + + Cl) (0-2 x + 0-2 2 + c 2 ) ’ 
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Method 

H <— constant 

A <— adaptive 

H <— adaptive 

A <— constant 

Figure 

Example 

Noise 

Level 

Stopping 

Parameter 

Convergence 

Time(s) 

PSNR (dB)/MSSIM 




14 (b) 

cr = 20 

e = 10- 4 

30.47 

28.7938/0.7723 

CTE (Proposed) 



14 (c) 

cr = 20 

e = 10 -6 

38.71 

27.3510/0.7147 




15 (b) 

cr = 30 

CD 

1 

O 

II 

38.61 

26.5238/0.6774 




14 (d) 

cr = 20 

e = 10- 4 

30.96 

28.2464/0.7237 

CTE (Proposed) 



14 (e) 

cr = 20 

e = 10 -6 

32.20 

26.1983/0.6796 




15 (c) 

cr = 30 

CD 

1 

O 

II 

32.33 

26.7451/0.6792 




14 (f) 

cr = 20 

e = 10- 4 

5.07 

27.3646/0.7178 

[48] 



14 (g) 

cr = 20 

e = 10 -6 

40.57 

25.4476/0.6970 




15 (d) 

cr = 30 

e = 10 -6 

40.30 

25.0559/0.6543 


Table 2: Image denoising error metrics (average) comparison using original coupled PDE scheme [48] for 
different noise levels and parameters on the Berkeley segmentation dataset (BSDS) 500. Some examples 
corresponding to the entries are shown above in Figure 14 and Figure 15. 

where /a UJi the average of c^, a^. the variance of u)i, cy u , lUJ2 the covariance, ci,C 2 stabilization pa¬ 
rameters, see [57] for more details. 

Remark 4. Note that the SSIM is a better error metric than PSNR as it provides a quantitative way 
of measuring the structural similarity of denoised image against the original noise-free image. To adapt 
PSNR and SSIM error metrics to color images one can convert the color image into gray 10 and then 
compute PSNR and MSSIM error metrics for the converted gray-scale image. In this paper, in order 
to compare with the scheme on [48], we compute the PSNR and MSSIM on each channel and use the 
average as a final value. 


4.3.2 Comparison with the previous model 

Table 2 compares proposed scheme with that of [48] using the both PSNR (dB) and MSSIM error metrics 
average for Berkeley segmentation dataset (BSDS) images. We implemented both the schemes on the 
full Berkeley segmentation dataset of 500 noisy images for two different noise levels and obtained similar 
improvements. Figures 14-15 show some example images corresponding to Table 2 . As can be seen from 
the zoomed in versions, the proposed CTE scheme provides cleaner cartoon components (denoised images, 
see Figures 14(b-e)-15(b,c)) in contrast to original coupled PDE model [48] which either excessively blurs 
out details (Figure 14(f,g)) or keeps noisy regions (Figure 15(d)) in final results. 

4.3.3 Comparison of adaptive fidelity parameters 

Figure 16 gives a comparison for the implementation of our CTE model (/i adaptive, A constant) 
by using adaptive fidelity parameters /ii (proposed), pL 2 (based on the work of [48]) and /i 3 (studied 
10 For example, using MATLAB’s command rgb2gray. 
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(a) Input noise (a = 20), (PSNR, MSSIM) = (24.1759,0.5130) 




(b) Proposed (ju <— adaptive, A -(—constant ) 


(c) Proposed (ju <— adaptive, A -(—constant) 


(PSNR, MSSIM) = (28.6484,0.7675) 


(PSNR, MSSIM) = (27.5146,0.6993) 




(d) Proposed (p, -(—constant, A -(—adaptive ) 


(e) Proposed (/x -(—constant, A -(—adaptive) 




(f) [48] (Adaptive) 

(PSNR, MSSIM) = (26.1172,0.7540) 


(g) [48] (Adaptive) 

(PSNR, MSSIM) = (24.8883,0.7057) 


Figure 14: (Color online) Better edge preserving image restoration results were obtained using our scheme 
in comparison with the original coupled PDE model [48] . The stopping parameters were chosen according 
to maximum MSSIM values, see Table 2 for the corresponding values. 
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(a) Input noise (a = 30), (PSNR, MSSIM) = (21.888,0.2445) 


(b) Proposed (p, <— adaptive, A constant) 


(PSNR, MSSIM) = (28.3610,0.7718) 




(c) Proposed (p, <— constant, A adaptive) 


(d) [48] (/x constant, A ^—adaptive) 


(PSNR, MSSIM) = (29.2489,0.7636) 


(PSNR, MSSIM) = (26.4529,0.6849) 


Figure 15: (Color online) Better edge preserving image restoration results were obtained using our scheme 
in comparison with the original coupled PDE model [48] . The stopping parameters were chosen according 
to maximum MSSIM values, see Table 2 for the corresponding values. 

in [30]) for a close-up flower image corrupted with Gaussian noise level a = 20. As shown different 
decomposition levels based on smooth + random noise + edges components are given according to the 
adaptive parameter /x. In terms of PSNR and MSSIM error metrics the fi i choice improves better the 
denoising result. Note the similar convergence property of the energy value between our /ii choice and 
the /i 3 adaptive parameter. 

4.3.4 Brain MRI image decomposition 

Figure 17 shows input Brain MR images from the Simulated Brain Database and its corresponding 
(u,v,w) functions for slices 50, 70 and 120 together with their corresponding inputs. Lack of sharp edges 
in the denoised Brain MRI images (Third-Fourth rows) can be attributed to the spatial relaxation based 
for the coupling PDE scheme. Spatial smoothing based regularization introduces a slight blur on each 
edge map w (Fifth row). Although, we gradually reduce the smoothing results by varying the stopping 
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(a) Input noise (cr = 20) 



(b) Energy vs Iteration for different /i choices 



(b) u,v,w,f — u components for our proposed scheme using \i i with (PSNR, MSSIM) = (27.6623, 0.8143) 


(c) u, v, w, f — u components for our proposed scheme using /Ji 2 with (PSNR, MSSIM) m (24.6084, 0.6823) 



(d) u, v, w, f — u components for our proposed scheme using ^,3 with (PSNR, MSSIM) = (25.6414, 0.7503) 



Figure 16: Comparison of different fi functions computed using a real image corrupted with Gaussian 
noise level a = 20 and stopping parameter e = 10 -4 . (a) Original image, (b) Energy versus iteration for 
different adaptive /i functions based energy minimization scheme (7-8). (c) fi 1 based on local histograms 
Eqn. (17). (d) /i 2 based on the work of [48] Eqn. (18). (e) /. 1 3 base on the work of [30] Eqn. (19). 
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(a) Noisy brain MRI images 


(c) Edges components (w) 


Figure 17: (Color online) u & w components using our proposed model with (/i adaptive, A constant) 
and stopping parameter e = 10 -4 for different Brain MR images (slice: 50, 70 and 120 respectively). 
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(b) Cartoon components (u) 







parameter e the strong discontinuities are preserved well and noise is removed effectively within regions. 


4.4 Multi-scale decomposition 

Following [54] we can make the weighted TV model with multi-scale parameter A. Note that this is 
slightly different from the original multi scale usage in [53], here we use it in the constraint Eqn. (8). Let 
us briefly recall the model proposed in [54] where the texture component v is modeled using G v norm. 
That is, the minimization is carried out for both u, v, 

^E(f,X;u,v) = \u\ B v(n) + Mil/- « - v\\ 2 L 2 (n) + A|M| Gp (n)j, 
with G p consisting of all distributions which can be written as, 

V = d x gi + d y g 2 = div(g), g € L P (Q,M. 2 ). 


The G p norm is defined as, 

IMIg p ( 0 ) = inf {11 <71/*>(£}) \v = div(g), geL p (fl,R 2 )} . 


We utilize the same modeling for the texture component v in our splitting step of the proposed weighted 
TV model (see Eqn. (4)), 


min 

u,v 


J ^ rVu| dx + Mil/ -u- vf L 2 (n) + A||v||g p (£ 


( 101 ) 


Finally, we compare multiscale version (101) with the multi-scale TV decomposition results of [54]. 
Figure 18 shows the comparison result on a synthetic image for 5 steps and our scheme retains the cartoon 
component clearly than the Tang and He [54]. Moreover, the texture components show a progressive 
capture of small scale oscillations. 


5 Conclusion 

We have presented a new image decomposition model coupling a variational and PDE via a weighted 
total variation regularization algorithm (CTE model). Our main contribution is twofold: 

1) The proposed decomposition model gets an image decomposition into its cartoon, texture and edge 
components with fixed, adaptive and multicale parameters for the id-fitting term by utilizing a local 
histogram along with a diffusion equation. Extensive experiments using a fast dual minimization 
based splitting implementation indicates that the proposed scheme is useful for edge preserving 
image decomposition on real and noisy images. Comparative results indicate the proposed scheme 
is useful in denoising natural images as well in multi scale image decomposition. 
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Multiscale cartoon and texture extraction using our proposed CTE method 



Multiscale cartoon and texture extraction using the hierarchical decomposition of Tang and He [54] 



Figure 18: (Color online) Multiscale cartoon and texture decomposition of a synthetic image for 5 steps by 
implementing our proposed approach (top two rows) (stooping parameter e = 10 -4 ) and the hierarchical 
method from [54] (bottom two rows). 
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2) We fashioned a new well posed scheme to transform the non-linear problem to a set of subproblems 
that are much easier to solve quickly via an optimization technique and a linear diffusion PDE 
solution. 

It is well known that there is no unique decomposition of an image into three scales: cartoon, texture 
and edges. At close range, texture and edges may be just a set of well-structured objects that can be 
kept in the cartoon part or the textural part according to a scale (equivalently iteration) parameter. We 
proposed an adaptive choice computed on local windows based histogram information as well as a multi¬ 
scale adaptation weight parameter for the id-fitting term, allowing us to get satisfactory results. There 
are other adaptive parameter choices for constructing a general decomposition model and we believe 
the proposed method is general in the sense that other regularizes (instead of TV) can be utilized in 
computing the cartoon component. We also remark that our denoising results are proof-of-concept for the 
proposed CTE model and we do not claim it outperforms state-of-the-art methods such as the nonlocal 
means [10] or BM3D [19] which are specifically designed for optimal denoising results. 
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